/* This file contains replication code for Vossler et al., "VALUING IMPROVEMENTS IN ECOLOGICAL INTEGRITY IN LOCAL AND REGIONAL WATERSHEDS" */
/* The associated data file is "BCG survey replication data.dta" */
/* For comments and questions, contact Christian Vossler (cvossler@utk.edu) */
/* April 2022 */

// Note: estimation relies on the user-written "mixlogit" command. You may have to install this command. 
// Note: the "outreg" command used to save results as an RTF file. You may have to install this command (or refrain from executing commands involving "outreg")
// Note: the "asdoc" command used to save results as a Word file. You may have to install this command (or refrain from executing commands involving "asdoc")

/* load data */
clear all
cd "D:\Research\Water Quality Valuation\Phase 1 analysis\"
use "Data\BCG survey replication data.dta"
/* initialize log file (optional) */
log using "Results\PNAS manuscript.log", replace

/* Table S4. Model 1, all voting scenarios */

// Estimate vacuous model; needed to compute R^2
clogit Vote, group(ChoiceID) cluster(RespondentID) 
scalar ll_0=e(ll)

// mixed logit in preference space; fixed cost parameter (Reported in Table S4)
mixlogit Vote Cost, rand(ASC WQ_HUC8 WQ_HUC4 WQ_Medium WQ_Large WQ_HUC4_NL WQ_Medium_NL) group(ChoiceID) id(RespondentID) cluster(RespondentID) nrep(500)
est store M1_FixedC_n500
// Calculate McFadden's R^2
display 1 - e(ll)/ll_0
scalar R2=1-e(ll)/ll_0
// Save model results
outreg2 using "Results\M1_FixedC_n500.rtf", stat(coef se) alpha (0.01, 0.05, 0.10) addstat(Log-likelihood, e(ll), McFadden's R2, R2) replace 

/* Table 2. Willingness-to-pay for selected water quality improvement scenarios */
// These results are derived from Model 1

// Calculate WTP for local sub-watershed (HUC8)
// Increase one BCG level across policy area
margins, express( -((_b[ASC] + _b[WQ_HUC8]*(HUC8BCG-1)) - (_b[WQ_HUC8]*HUC8BCG))/_b[Cost])
// Minimum BCG Level 2 ("swimmable")
margins, express( -((_b[ASC] + _b[WQ_HUC8]*2) - (_b[WQ_HUC8]*HUC8BCG))/_b[Cost])
// Minimum BCG Level 3 ("biological")
margins, express( -((_b[ASC] + _b[WQ_HUC8]*min(HUC8BCG, 3)) - (_b[WQ_HUC8]*HUC8BCG))/_b[Cost])

// Calculate WTP for local watershed (HUC4)
// Increase one BCG level across policy area
margins, express( -((_b[ASC] + _b[WQ_HUC4]*(HUC4BCG_base-1) + _b[WQ_HUC8]*(HUC8BCG-1)) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])
// Minimum BCG Level 2 ("swimmable")
margins, express( -((_b[ASC] + _b[WQ_HUC4]*HUC4BCG_scenario1 + _b[WQ_HUC8]*2) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])
// Minimum BCG Level 3 ("biological")
margins, express( -((_b[ASC] + _b[WQ_HUC4]*HUC4BCG_scenario2 + _b[WQ_HUC8]*min(HUC8BCG, 3)) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])

// Calculate WTP for local 3 Watersheds (3 HUC4s)
// Increase one BCG level across policy area
margins, express( -((_b[ASC] + _b[WQ_Medium]*(MediumBCG_base-1) + _b[WQ_HUC8]*(HUC8BCG-1)) - (_b[WQ_Medium]*MediumBCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])
// Minimum BCG Level 2 ("swimmable")
margins, express( -((_b[ASC] + _b[WQ_Medium]*MediumBCG_scenario1 + _b[WQ_HUC8]*2) - (_b[WQ_Medium]*MediumBCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])
// Minimum BCG Level 3 ("biological")
margins, express( -((_b[ASC] + _b[WQ_Medium]*MediumBCG_scenario2 + _b[WQ_HUC8]*min(HUC8BCG, 3)) - (_b[WQ_Medium]*MediumBCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])

// Calculate WTP for entire Study Region
// Increase one BCG level across policy area
margins, express( -((_b[ASC] + _b[WQ_Large]*(2.52) + _b[WQ_HUC8]*(HUC8BCG-1)) - (_b[WQ_Large]*3.52 + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])
// Minimum BCG Level 2 ("swimmable")
margins, express( -((_b[ASC] + _b[WQ_Large]*2 + _b[WQ_HUC8]*2) - (_b[WQ_Large]*3.52 + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])
// Minimum BCG Level 3 ("biological")
margins, express( -((_b[ASC] + _b[WQ_Large]*(3.52-0.59) + _b[WQ_HUC8]*min(HUC8BCG, 3)) - (_b[WQ_Large]*3.52 + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])

// Calculate WTP for non-local watershed (HUC4)
// Increase one BCG level across policy area
margins, express( -((_b[ASC] + _b[WQ_HUC4_NL]*(HUC4BCG_base_NL-1) ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost])
// Minimum BCG Level 2 ("swimmable")
margins, express( -((_b[ASC] + _b[WQ_HUC4_NL]*2 ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost])
// Minimum BCG Level 3 ("biological")
margins, express( -((_b[ASC] + _b[WQ_HUC4_NL]*HUC4BCG_scenario2_NL ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost])

// Calculate WTP for non-local 3 Watersheds (3 HUC4s)
// Increase one BCG level across policy area
margins, express( -((_b[ASC] + _b[WQ_Medium_NL]*(MediumBCG_base_NL-1) ) - (_b[WQ_Medium_NL]*MediumBCG_base_NL ))/_b[Cost])
// Minimum BCG Level 2 ("swimmable")
margins, express( -((_b[ASC] + _b[WQ_Medium_NL]*2 ) - (_b[WQ_Medium_NL]*MediumBCG_base_NL ))/_b[Cost])
// Minimum BCG Level 3 ("biological")
margins, express( -((_b[ASC] + _b[WQ_Medium_NL]*MediumBCG_scenario2_NL ) - (_b[WQ_Medium_NL]*MediumBCG_base_NL))/_b[Cost])

/* Some hypothesis tests based on Model 1 */
// HUC4 v. Medium local, min2 scenario
margins, express(( -((_b[ASC] + _b[WQ_HUC4]*HUC4BCG_scenario1 + _b[WQ_HUC8]*2) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost]) - ( -((_b[ASC] + _b[WQ_Medium]*MediumBCG_scenario1 + _b[WQ_HUC8]*2) - (_b[WQ_Medium]*MediumBCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost]))
// HUC4 v. Large local, min2 scenario
margins, express(( -((_b[ASC] + _b[WQ_HUC4]*HUC4BCG_scenario1 + _b[WQ_HUC8]*2) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost]) - ( -((_b[ASC] + _b[WQ_Large]*2 + _b[WQ_HUC8]*2) - (_b[WQ_Large]*3.52 + _b[WQ_HUC8]*HUC8BCG))/_b[Cost]))
// HUC4 v. Large local, 1-up scenario
margins, express(( -((_b[ASC] + _b[WQ_HUC4]*(HUC4BCG_base-1) + _b[WQ_HUC8]*(HUC8BCG-1)) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost]) - ( -((_b[ASC] + _b[WQ_Large]*(2.52) + _b[WQ_HUC8]*(HUC8BCG-1)) - (_b[WQ_Large]*3.52 + _b[WQ_HUC8]*HUC8BCG))/_b[Cost]))

// HUC4 v. Medium non-local, min2 scenario
margins, express( (-((_b[ASC] + _b[WQ_HUC4_NL]*2 ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost]) - ( -((_b[ASC] + _b[WQ_Medium_NL]*2 ) - (_b[WQ_Medium_NL]*MediumBCG_base_NL ))/_b[Cost]))
// HUC4 v. Medium non-local, min3 scenario
margins, express( (-((_b[ASC] + _b[WQ_HUC4_NL]*HUC4BCG_scenario2_NL ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost]) - ( -((_b[ASC] + _b[WQ_Medium_NL]*MediumBCG_scenario2_NL ) - (_b[WQ_Medium_NL]*MediumBCG_base_NL ))/_b[Cost]))
// HUC4 v. Medium non-local, 1-up scenario
margins, express(( -((_b[ASC] + _b[WQ_HUC4_NL]*(HUC4BCG_base_NL-1) ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost]) - ( -((_b[ASC] + _b[WQ_Medium_NL]*(MediumBCG_base_NL-1) ) - (_b[WQ_Medium_NL]*MediumBCG_base_NL ))/_b[Cost]))

/* Some robustness checks related to Model 1 */

// Robustness check: conditional logit
clogit Vote ASC WQ_HUC8 WQ_HUC4 WQ_Medium WQ_Large WQ_HUC4_NL WQ_Medium_NL Cost, group(ChoiceID) cluster(RespondentID) 
est store M1_clogit
// Calculate McFadden's R^2
display 1 - e(ll)/ll_0
scalar R2=1-e(ll)/ll_0

// Robustness check: mixed logit with sampling weights
mixlogit Vote Cost [pweight=WEIGHT], rand(ASC WQ_HUC8 WQ_HUC4 WQ_Medium WQ_Large WQ_HUC4_NL WQ_Medium_NL) group(ChoiceID) id(RespondentID) cluster(RespondentID) nrep(500)
est store M1_FixedC_n500_wt
// Calculate McFadden's R^2
display 1 - e(ll)/ll_0
scalar R2=1-e(ll)/ll_0

// Robustness check: mixed logit with correlated parameters
mixlogit Vote Cost, rand(ASC WQ_HUC8 WQ_HUC4 WQ_Medium WQ_Large WQ_HUC4_NL WQ_Medium_NL) group(ChoiceID) id(RespondentID) cluster(RespondentID) nrep(500) corr
est store M1_FixedC_n500_corr
// Calculate McFadden's R^2
display 1 - e(ll)/ll_0
scalar R2=1-e(ll)/ll_0


/* Table S5. Model 2, local and non-local watershed voting scenarios */

// Estimate vacuous model; needed to compute R^2
clogit Vote if (Spatial_Unit==1 | Spatial_Unit==4), group(ChoiceID) cluster(RespondentID) 
scalar ll_0=e(ll)

// mixed logit in preference space; fixed cost parameter
mixlogit Vote ASC_PercentInstate ASC_PercentInstate_NL Cost if (Spatial_Unit==1 | Spatial_Unit==4), rand(ASC WQ_HUC8 WQ_HUC4 WQ_HUC4_NL) group(ChoiceID) id(RespondentID) cluster(RespondentID) nrep(500)
est store M2_FixedC_n500
// Calculate McFadden's R^2
display 1 - e(ll)/ll_0
scalar R2=1-e(ll)/ll_0
// Save model results
outreg2 using "Results\M2_FixedC_n500.rtf", stat(coef se) alpha (0.01, 0.05, 0.10) addstat(Log-likelihood, e(ll), McFadden's R2, R2) replace 


/* Table 3. Willingness-to-pay for water quality improvement scenarios based on percentage of impacted area located in-state */
// These results are derived from Model 2

// Local policy: impact area 100% in-state
margins, express( -((_b[ASC] + _b[ASC_PercentInstate]*100 + _b[WQ_HUC4]*(HUC4BCG_base-1) + _b[WQ_HUC8]*(HUC8BCG-1)) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])
margins, express( -((_b[ASC] + _b[ASC_PercentInstate]*100 + _b[WQ_HUC4]*HUC4BCG_scenario1 + _b[WQ_HUC8]*2) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])
margins, express( -((_b[ASC] + _b[ASC_PercentInstate]*100 + _b[WQ_HUC4]*HUC4BCG_scenario2 + _b[WQ_HUC8]*min(HUC8BCG, 3)) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])

// Local policy: impact area 25% in-state
margins, express( -((_b[ASC] + _b[ASC_PercentInstate]*25 + _b[WQ_HUC4]*(HUC4BCG_base-1) + _b[WQ_HUC8]*(HUC8BCG-1)) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])
margins, express( -((_b[ASC] + _b[ASC_PercentInstate]*25 + _b[WQ_HUC4]*HUC4BCG_scenario1 + _b[WQ_HUC8]*2) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])
margins, express( -((_b[ASC] + _b[ASC_PercentInstate]*25 + _b[WQ_HUC4]*HUC4BCG_scenario2 + _b[WQ_HUC8]*min(HUC8BCG, 3)) - (_b[WQ_HUC4]*HUC4BCG_base + _b[WQ_HUC8]*HUC8BCG))/_b[Cost])

// Non-local policy: impact area 25% in-state
margins, express( -((_b[ASC] + _b[ASC_PercentInstate_NL]*25+ _b[WQ_HUC4_NL]*(HUC4BCG_base_NL-1) ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost])
margins, express( -((_b[ASC] + _b[ASC_PercentInstate_NL]*25 + _b[WQ_HUC4_NL]*2 ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost])
margins, express( -((_b[ASC] + _b[ASC_PercentInstate_NL]*25+ _b[WQ_HUC4_NL]*HUC4BCG_scenario2_NL ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost])

// Non-local policy: impact area 0% in-state
margins, express( -((_b[ASC]  + _b[ASC_PercentInstate_NL]*0 + _b[WQ_HUC4_NL]*(HUC4BCG_base_NL-1) ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost])
margins, express( -((_b[ASC]  + _b[ASC_PercentInstate_NL]*0 + _b[WQ_HUC4_NL]*2 ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost])
margins, express( -((_b[ASC]  + _b[ASC_PercentInstate_NL]*0 + _b[WQ_HUC4_NL]*HUC4BCG_scenario2_NL ) - (_b[WQ_HUC4_NL]*HUC4BCG_base_NL ))/_b[Cost])

/* Some robustness checks related to Model 2 */

// Robustness check: conditional logit
clogit Vote ASC_PercentInstate ASC_PercentInstate_NL WQ_HUC8 WQ_HUC4 WQ_HUC4_NL Cost if(Spatial_Unit==1 | Spatial_Unit==4), group(ChoiceID) cluster(RespondentID) 
est store M2_clogit
// Calculate McFadden's R^2
display 1 - e(ll)/ll_0
scalar R2=1-e(ll)/ll_0

// Robustness check: mixed logit with sampling weights
mixlogit Vote ASC_PercentInstate ASC_PercentInstate_NL Cost [pweight=WEIGHT] if (Spatial_Unit==1 | Spatial_Unit==4), rand(ASC WQ_HUC8 WQ_HUC4 WQ_HUC4_NL) group(ChoiceID) id(RespondentID) cluster(RespondentID) nrep(500)
est store M2_FixedC_n500_wt
// Calculate McFadden's R^2
display 1 - e(ll)/ll_0
scalar R2=1-e(ll)/ll_0

// Robustness check: mixed logit with correlated parameters
mixlogit Vote ASC_PercentInstate ASC_PercentInstate_NL Cost if (Spatial_Unit==1 | Spatial_Unit==4), rand(ASC WQ_HUC8 WQ_HUC4 WQ_HUC4_NL) group(ChoiceID) id(RespondentID) cluster(RespondentID) nrep(500) corr
est store M2_FixedC_n500_corr
// Calculate McFadden's R^2
display 1 - e(ll)/ll_0
scalar R2=1-e(ll)/ll_0


/* Table S6. Descriptive statistics for selected survey questions and socio-economic characteristics */

sum i.P3_Voted_WouldPay i.P3_Voted_OutcomeAchieved i.P3_Voted_InformPolicy i.P3_Importance_Size i.P3_Importance_WQ i.P3_Importance_Cost P4_Recreation_Trips i.P4_FarthestTrip_Cat ///
Q_Female Q_Age i.Q_Race i.Q_Educ_Cat Q_Married Q_Retired Q_Income Q_Metro Q_HHsize if Referendum==1 & ASC==1

asdoc sum i.P3_Voted_WouldPay i.P3_Voted_OutcomeAchieved i.P3_Voted_InformPolicy i.P3_Importance_Size i.P3_Importance_WQ i.P3_Importance_Cost P4_Recreation_Trips i.P4_FarthestTrip_Cat ///
Q_Female Q_Age i.Q_Race i.Q_Educ_Cat Q_Married Q_Retired Q_Income Q_Metro Q_HHsize if Referendum==1 & ASC==1, replace save(Results\Summary.doc) title(Descriptive Statistics) 


/* Table S7. Effects of socio-economic characteristics on the willingness-to-pay for water quality improvements */ 

global ASC_covariates ASC_Q_Female ASC_Q_Age ASC_Q_Educ_2 ASC_Q_Educ_3 ASC_Q_Educ_4 ASC_Q_Educ_5 ASC_Q_Race_2 ASC_Q_Race_3 ASC_Q_Race_4 ASC_Q_Race_5 ASC_Q_Race_6 ASC_Q_Married ASC_Q_Retired ASC_Q_Income ASC_Q_Metro ASC_Q_HHsize

// Estimate vacuous model; needed to compute R^2
clogit Vote, group(ChoiceID) cluster(RespondentID) 
scalar ll_0=e(ll)

// mixed logit in preference space; fixed cost parameter
mixlogit Vote Cost $ASC_covariates, rand(ASC WQ_HUC8 WQ_HUC4 WQ_Medium WQ_Large WQ_HUC4_NL WQ_Medium_NL) group(ChoiceID) id(RespondentID) cluster(RespondentID) nrep(500)
est store M3_FixedC_Demo
// Calculate McFadden's R^2
display 1 - e(ll)/ll_0
scalar R2=1-e(ll)/ll_0
// Save model results
outreg2 using "Results\M3_FixedC_Demo.rtf", stat(coef se) alpha (0.01, 0.05, 0.10) addstat(Log-likelihood, e(ll), McFadden's R2, R2) replace 

// convert utility parameters into willingness-to-pay estimates (as reported in table)
foreach var in $ASC_covariates {
nlcom -_b[`var']/_b[Cost]
}

// test for equality of race/ethnicity categories
test (ASC_Q_Race_2 =ASC_Q_Race_3 =ASC_Q_Race_4 =ASC_Q_Race_5 =ASC_Q_Race_6=0)

log close


